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Abstract We study properties of chaos in generic one-dimensional nonlinear Hamiltonian lattices 
comprised of weakly coupled nonlinear oscillators, by numerical simulations of continuous-time systems 
and symplectic maps. For small coupling, the measure of chaos is found to be proportional to the 
coupling strength and lattice length, with the typical maximal Lyapunov exponent being proportional 
to the square root of coupling. This strong chaos appears as a result of triplet resonances between 
nearby modes. In addition to strong chaos we observe a weakly chaotic component having much smaller 
Lyapunov exponent, the measure of which drops approximately as a square of the coupling strength 
down to smallest couplings we were able to reach. We argue that this weak chaos is linked to the 
regime of fast Arnold diffusion discussed by Chirikov and Vecheslavov. In disordered lattices of large 
size we find a subdifiusive spreading of initially localized wave packets over larger and larger number 
of modes. The relations between the exponent of this spreading and the exponent in the dependence of 
the fast Arnold diffusion on coupling strength are analyzed. We also trace parallels between the slow 
£SJ ' spreading of chaos and deterministic rheology. 
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1 Introduction 



Even 120 years after the fundamental work of Poincare [T] and numerous efforts done after it, an 
interplay between order and chaos in high-dimensional Hamiltonian systems remains a challenging 
problem. For Hamiltonian dynamics with a few degrees of freedom, a clear picture of a separation 
between chaotic and regular (quasiperiodic) regions in the phase space [2l[3] has been confirmed in 
numerous studies. Much less is known on this separation and the structural properties of chaos when 
the number of degrees of freedom becomes large. Especially the generic case of a weak nonlinear 
coupling of initially nonlinear but integrable degrees of freedom remains poorly understood. We will 
call such systems to be weakly nonintegrable. Their properties are very nontrivial since a decrease 
in nonlinearity/nonintegrability might be compensated by an increase of dimensionality of the phase 
space. 

The Kolmogorov-Arnold-Moser (KAM) theory guarantees the existence of invariant tori at a suffi- 
ciently weak nonlinear perturbation (see e.g. [2113] )■ However, in conservative systems with more than 
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two degrees of freedom (L > 2) such tori are not isolating and chaos can spreads along tiny chaotic 
layers as it was pointed by Arnold [3]. The mechanism of such a chaotic spreading is known under 
the name of Arnold diffusion as coined by Chirikov in 1969 [2,5 . For L = 3 the rate of Arnold dif- 
fusion drops exponentially with the dimensionless strength of nonlinear coupling ft 2;S T 5 i . This is in 
a qualitative agreement with a number of mathematical results which give rigorous bounds on the 
spreading rate in the limit of asymptotically small ft at fixed L > 2 [5J[7]. The mathematical studies 
of the Arnold diffusion properties are actively continued at present (see e.g. [8] and Refs. therein). 
While the mathematical results indicate the exponentially small rate of Arnold diffusion Da in the 
limit of small nonlinearity ft at fixed L, it remains not clear at what realistic values of nonlinearity 
such an exponential behavior effectively appears. The striking results of Chirikov and Vecheslavov, 
established by extensive numerical simulations for 4 < L < 15 and supported by heuristic arguments 
9, lOlfTTj. show only an algebraic decay of Da with ft up to extremely small values of Arnold diffusion 
coefficient Da ~ 10 -50 . This regime was named by them the fast Arnold diffusion. These studies have 
been restricted by L < 15 and it remains unclear what can happen with such a behavior in the limit 
of larger L with small but fixed ft. 

The question about the properties of Hamiltonian systems at large values of L is linked to the fun- 
damental problem of dynamical thermalization and ergodicity in the thermodynamic limit. As typical 
models with a large number of degrees of freedom one considers Hamiltonian lattices (or Hamiltonian 
partial differential equations, which, however, live in an infinite-dimensional phase space). A striking 
example of nontrivial dynamics in weakly nonlinear lattices gives the Fermi-Pasta- Ulam problem |12j , 
which is still far from being completely resolved despite of numerous efforts in its 50-year history (see 
[T3lfl4] for a stand around 2004 and [15] for recent advances). Moreover, because the FPU model has a 
special peculiarity as being close to an integrable Toda lattice, its properties appear to be rather spe- 
cial. Quite recently, a lot of attention attracted disordered nonlinear lattices [16, 17, 18, 19,20,21 ,22,23, 
24,25,26,27 studied in the context of the problem of nonlinear destruction of Anderson localization. 
Here one tries to relate the properties of chaos and regularity at small nonlinearities to the properties 
of the spreading of a wave packet over the lattice [2"51l2"9"lf3"0"] . Certain mathematical bounds on the rate 
of spreading have been obtained [31 , 32 by the methods similar to those of Nekhoroshev [B] but they 
are available only in the limit of very small nonlinearity being very far from the regimes studied in 
numerical simulations. In addition, these weakly nonlinear lattices are not generic objects from the 
point of view of weak nonintegrability and the KAM theory, since in the limit of small coupling they 
are reduced to a set of linear modes, i.e. to a linear quasiperiodic state demonstrating qusiperiodicity 
and pure point spectrum typical of the Anderson localization, and not to the generic case with a set 
of uncoupled nonlinear modes. We note, that in context of the KAM theory, a small perturbation of 
the latter integrable nonlinear system is studied. 

In this paper we study properties of a lattice of weakly coupled nonlinear oscillators at small cou- 
pling and large number of degrees of freedom. In the limit of small coupling this model reduces to an 
integrable although strongly nonlinear one, demonstrating typically quasiperiodic dynamics. A nice 
model of such a setup has been suggested by Kaneko and Konishi [331134] . it gives a generalization of 
the Chirikov standard map [2] to a lattice of coupled symplectic maps. This model is computationally 
efficient and allows one a rather good numerical characterization of properties of regular and chaotic 
dynamics. Nevertheless, even for this model the quantitative properties are not well-established de- 
spite of various efforts f35"]13BHT0" |l 1 1 [157] . Additionally, we study here two models of coupled nonlinear 
continuous-time oscillator lattice where the spreading over the lattice can be analyzed at fixed energy. 
Our main conclusions are valid for all these systems. 

The plan of the paper is as follows. We start by formulating basic models we study in Section [5] 
Then in Section [3] we discuss the properties of the largest Lyapunov exponent, especially the scaling 
relations in dependence on coupling strength and system length. In Section U we argue that chaos is 
mainly due to occasional resonances between triples of three neighboring oscillators. In Section [5] we 
discuss statistical properties of chaos, focusing on the scaling of the diffusion constant. In Section H2 
we relate this properties to that of spreading of a wave packet in an unbounded lattice. Finally, in 
Section [7J a very slow evolution is compared to similar effects in the context of rheology. 
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2 Basic models 

Here we introduce three generic models of nonlinear oscillators locally coupled in space. Model A, 
introduced by Kaneko and Konishi [33,34 ], is a model of coupled symplectic maps 

Pk = Pk + K[sm(xk+i - xu) + sin(x fe _i - x k )] , k = l,...,L 

X k = Xk+Pk ■ 

with periodic boundary conditions. Here p is a "momentum" and a; is a "phase" variable. In the 
absence of coupling (i.e. for K = 0) each oscillator has a constant frequency pk that depends on initial 
conditions, so in the whole lattice generally a quasiperiodic regime with L frequencies establishes. For 
finite K the oscillators are coupled and chaos is possible. 

Model B is a strongly nonlinear continuous-time lattice with Hamiltonian 

# = + | + f(* +1 -«Z fe ) 2 . (2) 

fc=i 

Here we also consider a lattice of length L with periodic boundary conditions. The coupling parameter 
j3 plays the same role as K . Contrary to model A, model B conserves the total energy. We normalize 
the energy in such a way that E = L (i.e. density of energy is one), so that (3 and L are the only 
parameters of this model. 

Very similar to the model B is the model C, where the coupling between nonlinear modes is also 
nonlinear, moreover, the power of nonlinearity in coupling is stronger than the local one: 

H = t P j + ^ + ^ fe+1 -*) 6 , (3) 

k=l 

where we consider two cases for coefficients with all r\k — 1 (CI) and random homogeneously distributed 
values 0.5 < r\k < 1-5 (C2). While we do not expect large difference between models B and C in the 
described setup, where the density of the energy is fixed, the situation changes when the total energy 
is fixed and the length of the lattice is increased. In this limit model B will become asymptotically 
linear (effective j3 increases) while model C will become asymptotically less and less coupled (effective 7 
decreases). This difference is important for the implications of chaos for spreading of initially localized 
wave packets, to be discussed in Section [S] The randomness of values of r\k (model C2) ensures that 
there are no regular waves emanating from the main part of the wave packet in contrast to the case 
rjk = 1 (model CI) where such wave radiation is possible 39, 4"0]. 



3 Lyapunov exponents and their scaling 

3.1 Lyapunov exponents 

The largest Lyapunov exponent (LE) is a standard measure of chaos and is easy to calculate [31155] . We 
have performed a statistical analysis of Lyapunov exponents for models A, B, C based on an ensemble 
of random initial conditions. For model A we have chosen < pk,Xk < 27r as independent uniformly 
distributed. For model B we initialized qk = and pk normally distributed with zero mean, after this 
the values pk are rescaled such that the total energy of the lattice equals L - the number of lattice 
sites. For the model C the initialization is done in a similar way. We used up to several thousands of 
initial state realizations to obtain a good statistics in the computation of measure of chaos P c h- 

We present the "raw data" of these calculations for models A and B in Fig. [T] Here, for model A 
in a lattice with L = 8 one observes predominantly chaos for K = 0.05, predominantly regularity for 
K = 0.001, and both states depending on initial conditions for K = 0.01. Noteworthy, LE in the case 
of regularity does not vanish but attains very small values, with the cutoff appearing due to a finite 
integration time. In the middle part of Fig. HJa) one can see that increasing the integration time by 
factor 10 roughly decreases this lower cutoff in the Lyapunov exponent by factor 10. For any fixed T av , 
basing on inspection, one easily chooses a threshold in LE that separates chaos from regularity. Of 
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course, there are realizations with values around these thresholds that cannot be resolved within the 
integration time used, but their statistical relevance is not significant. Essentially the same picture is 
observed for models B (Fig. [TJd) and model C (data not shown). 
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Fig. 1 (Color online) (a) Calculations of LEs for model A with L = 8. First 10000 points (red): K = 0.001, LE 
calculated by averaging over time interval T av — 5 • 10 6 . Second 10000 points (green): K = 0.01, T av — 5 • 10 6 . 
Next 10000 points (also green): K = 0.01, T av = 5 ■ 10 7 . Last 10000 points (blue): K = 0.05, T av = 5 • 10 6 . 
(b): The same as (a) but for model B with L — 16, T av — 10 6 and different /3, from left to right: /3 = 
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for each value of (3 3000 realizations are shown. 



3.2 Scaling of probability to observe chaos 

According to calculations of LEs we can distinguish chaotic and regular regimes, and calculate the 
probability of their appearance in models A, B, C. The results for coupled symplectic maps of model 
A are presented in Fig. [5J A typical lower cutoff for the LE calculated over time interval T = 10 8 was 
« 2.5 • 10~ 8 , so we attributed all the realizations with A > 5 ■ 10~ 7 to chaos. Defined in this way the 
total measure of initial conditions in the phase space that yield chaos P c h decreases with K and L. 
The rescaled plot shows that for small K and large L the scaling relation 

Pch-K-L (4) 

holds. The same scaling P c h ~ K ■ f3 is valid also for model B, as demonstrated in Fig. El an d for model 

c (Fig.®. 

The scaling with the system length P c h ~ L has been already discussed for model A in [351136] and for 
disordered nonlinear lattices in [3D] - It is based on the locality of chaos: the latter appears due to a local 
in space nonlinear interaction of localized modes, and not due to propagation of waves. Thus, in order 
to observe regularity in the whole lattice, the dynamics has to be regular in all subparts. Therefore, if 
the measure of chaos in a sublattice of length L , P c h{Lo), is small, then P c h(L) « 1 — (1 — P c h(Lo)) L ^ L ° 
from which the scaling log P c h ~ L follows. An additional check of this relation is in Fig. [SJd below. 
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Fig. 2 Calculations of Pch vs K for model A, demonstrating the scaling P c h ~ KL for small K. The middle 
panel shows the same data as the left one but in a logarithmic scale, while the left panel shows P c h as a function 
of the product KL. The dashed line on the right panel is P c h = KL. 
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Fig. 3 The same as Fig. [2] but for model B. 




Fig. 4 The same as Fig. 0but for model CI. 
3.3 Scaling of the Lyapunov exponent 

Next, we studied the scaling properties of the value of LE. As one can already see from Fig. [TJ the 
positive LEs concentrate around a maximal value that decreases with K and /3. We have found (see 
Fig. [5^,) that this maximal value is roughly independent on the length of the system L and scales with 
nonlincarity parameters K and j3 as 

A ~ K 1 ' 2 . (5) 




Fig. 5 Left panel: Dependence of the Lyapunov exponent A on the perturbation parameter K for model A 
at L = 4 (blue circles), 8 (red squares). Dashed straight line shows approximate dependence A = 0.1Q\/~K (for 
N = 4, 8). The fit of data gives the exponent of the dependence A oc K a with a = 0.476 ± 0.013 (for L = 4), 
0.450 ± 0.012 (for L = 8) in agreement with the scaling |5|. Right panel: Dependence of the measure of chaos 
P c h on the perturbation parameter K for L — 4 (blue circles), 8 (red squares). Dashed straight line shows 
approximate dependence P c h ~ 7.75K (for L — 4); for L — 8 we find that P c h = 14.32K in agreement with the 
scaling Q. The fit of data gives the exponent of the dependence P c h oc K b with b — 1.022 ± 0.01 (for L = 4), 
1.036 ± 0.01 (for L = 8). Up to 5 x 10 5 trajectories and time t < 10 6 have been used to compute the Lyapunov 
exponent A and determine the number of chaotic trajectories with A > 0. Certain checks have been made with 
t — 5 x 10 9 and 100 trajectories. 

To demonstrate the scaling of the Lyapunov exponents we calculated their probability distribution 
densities w(X). Because of the relation P c h — iv(X) dX (where X t h is the cutoff value) the appro- 
priate scaling for this density is that of P c h, i.e. K ■ L. According to ©, the appropriate scaling of 
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the argument of the density is XK~ X I 2 . We plot rescaled in this way distribution densities of LEs for 
models A and B in Figs. H2 and [7] We present here results for the distribution density w, for constructing 
of which some arbitrary bins have been used, and for a cumulative distribution W(A) = f. w(\) dX 
where all data are presented, respectively. We note that the scaling law differs from the scaling 
A ~ K 2 / 3 suggested in [351136] . For the model B we find the same scaling relation A ~ /3 1 / 2 as it is 
shown in Fig. [7}d. For the model C we find the similar relation. 
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(a) log^A-tf- 1 / 2 ) ( b ) log 10 (A-r 1/2 ) 

Fig. 6 (a): Distribution density of LEs in model A. (b): Distribution density of LEs in model B. 
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Fig. 7 (a): cumulative distribution of LEs in model A. (b): cumulative distribution of LEs in model B. 



3.4 Strong and weak chaos 

There is also a substantial part of trajectories that have LEs between the lowest cutoff (determined 
by the averaging time) and the largest value ~ \K. We will distinguish these regimes by referring to 
the dynamics with LEs in the peak of distribution in Fig. [6] as strong chaos while the dynamics with 
lower LEs will be called weak chaos. As it will be discussed later, it might be that the regime of weak 
chaos is that where the fast Arnold diffusion discussed in [TT] occurs. We show in Fig. [5] that the total 
probability P wc h to observe this weak chaos scales as 

P wch oc K"™ ch L , v wch w 1.6 . (6) 

In Fig. IH1 we show an example of a local in time LEs for one long trajectory in model A. It shows 
existence of transitions between regimes with strong chaos and weak chaos. 
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K 

Fig. 8 Probability of "weak chaos" P wch with the LE between the low A > 10~ 7 and high (A < 0.03 ■ K 1/2 ) 
cutoffs in Fig.[7|in model A. The pink dotted line is P wc h ~ K 1 ' 6 , the dashed curve corresponds to the estimate 
Pwch ~ K 2 5 (lnK) 2 , discussed at the end of Section[S]in relation to the regime of fast Arnold diffusion analyzed 
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Fig. 9 Local LEs in model A as a function of time. Each value of LE is calculated over time interval of length 
10 6 , here 2 • 10 3 such time intervals are shown; K — 0.05, L = 8. 



4 Resonances as a source of chaos 

In order to characterize conditions under which chaos occurs at very small coupling, we have looked 
on resonances, and have found that chaos is highly correlated with the triple resonance at which the 
frequencies of three neighboring oscillators nearly coincide. For models A and B we illustrate this in 
Figs. [TUJ respectively. Here the probability of chaos P c h is shown vs. renormalized distances of initial 
frequencies of oscillators. For model A we have defined this distance as d = rrnnfc[(/(pjj. — pjK^) 2 + 
(/(Pfc+i — Pfc+2)) 2 ]- Here f(x) — 2| sin0.5x| measures the closeness of two initial momenta modulo 2ir. 
A small value of d indicates that somewhere in the lattice three initial nearby momenta p^ are close 
to each other. Then, for different realizations of initial conditions, different K in the range [0.001, 0.2] 
and different lattice lengths L = 8, 16,32 we determined the probability for chaos to occur vs. d/\K. 
One can see that for different lattice lengths the curves are close to each other, thus indicating that 
indeed the occurrence of resona nces is a nec essary prereq uisit e for cha os. In a similar analysis for model 
B we used d 2 = mm k [(y/p k (0) - y/ Pk+1 (0)) 2 + (y/p k (0) - v'Pfc-i(O)) 2 ]. 

In Fig. [TU] we demonstrate the correlation between the occurrence of resonance (small d) and the 
probability to observe chaos P c h- Moreover, we see here the scaling that in fact d should be compared 
with \[K (or \/]3 for model B). 

The physical reason for the scaling results presented in previous sections is the following (for 
simplicity of presentation, we refer here to model A only, the same arguments work for models B and 
C). There is a finite probability that three nearby particles will have their frequencies Ui — pi close to 
each other, within the frequency range Auj = \J~K. The probability of such an event is P ~ K, since 
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(a) d/VK (b) d/VP 

Fig. 10 The probability to observe chaos in dependence on the resonance in initial data for models A (a) and 
B (b). The data for different L and K,/3 collapse if the distance in initial frequencies is scaled according to 
VKot 



the first particle may have any frequency, the probability to have the second in the range \[~K is \J~K 
and the probability to have the third in the same range is also \f~K. This gives the probability of the 
resonance P ~ K for a lattice with three particles and P ~ K L for a chain with L oscillators. Similar 
arguments work for models B,C. It is important to note that in the case of such a 3-particle resonance, 
the KAM arguments are not valid and the dynamics remains chaotic at arbitrary small perturbation 
K. The situation is similar to the one considered in 0T] where three linear oscillators with the same 
frequency remain chaotic at arbitrary small nonlinear coupling between them. Indeed, in our case the 
numerical analysis shows that almost all chaotic trajectories (those with positive Lyapunov exponent) 
have three nearby particles with close frequencies. 




-3-2-10123 -3-2-10123 

(a) (b) 4>2 

Fig. 11 Poincare sections of variables 02, Ji for the resonance Hamiltonian (|12|l at 0i — 0. (a): Hr = 0, here 
chaos is dominant, (b): Hr = 10, here the dynamics is typically quasiperiodic. 



To understand this phenomenon in a better way let us consider the case when initially at three 
neighboring sites the values of actions pi are close to their average value P = (pi + p 2 +ps)/3. Then 
the evolution of these three particles, considered separately from the rest (what can be justified by 
arguing that nonresonant terms effectively disappear after averaging) is described by the mapping 

pi -pi = Ksm(x 2 -xi) , xi-xi=pi, (7) 
P2 - P2 = Ksm(xi - x 2 ) + K sin(a; 3 - x 2 ) , x 2 - x 2 — p 2 , (8) 
P3 - P3 = K sm(x 2 - xa) , x 3 - x 3 = p 3 . (9) 

Exploring the integral p\ = const and performing a canonical transformation to new conjugate 

coordinates according to 

01 = Xl-X 2 , 02 = x 3 ~x 2 , 03 = xi+x 2 +x 3 , pi = Ii+I 2 + P, p 2 = -h-I 2 +I 3 + P, p 3 = I 2 +I 3 +P, 

we obtain a two-dimensional mapping 

h-h = -K sin 0i , 0i - 0i = 2Ji + I 2 , (10) 
h~h = ~K Sin 02 , 02 - ^2 = 2l 2 + h , (11) 
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which due to smallness of I\ , I2 and of K can be approximated as a continuous-time system with 
Hamiltonian H = I 2 + /f + — K cos^i — K cos fa. After rescaling of actions to Ji = U/^/K and 
time to t = \ Kt we come to dimensionless resonance Hamiltonian 

E R (Ji,J 2 , fa, fa) = J\ + J| + J1J2 - cos(fa) - cos(fa) (12) 

Note that this rescaling proofs the dependencies ~ K 1 / 2 for the allowed deviations from the resonance 
condition. Also the rescaling of time proofs the scaling of the Lyapunov exponent with K according 
to©. 

According to the Chirikov resonance-overlap criterion [2] the dimensionless dynamics of Hamilto- 
nian (fT2"|) is chaotic for small values of energy (i.e. close to resonance) and chaos disappears if the 
energy is large. The Poincare sections for Hr for Hr = and Hr = 10 are shown in Fig. [TT1 confirming 
this picture. 

5 Properties of diffusion and weak chaos 

While LEs serve as an important indication for chaos, other quantities like correlations are important 
to characterize irregularity of the dynamics. For the Chirikov standard map an important statistical 
quantity is the diffusion constant of the momentum p: at large times T the dynamics of p can be 
considered as a random walk with a diffusion constant D defined according to ((p(T) — p(0)) 2 ) = DT. 
For the Chrikov standard map the dependence of D on the parameter K is known in detail [HE]. 

For the coupled symplectic maps (model A) numerical computations 33,34 , performed in a range 
0.1 < K < 1, indicated a weak diffusion at K — 0.1, the authors fitted the data with a stretched 
exponential dependence. Here we extend these calculations and show the results in Fig. [12] One can 
see a strong decrease of the diffusion constant with K, which for small K is close to a power-law 
dependence 

D ~ K VD , v D w 6.5 . (13) 

A similar value of the exponent was obtained from the statistics of Poincare recurrences in the range 
0.1 < K < 1 [42]. We note that for model C the above equation implies D oc 7"°. The value of the 
exponent vr, is close to the value given by Chirikov and Vecheslavov [TUl[TTj . However, they calculated 
the diffusion indirectly by expressing it via an effective width w s of a separatrix layer of a nonlinear 
resonance with the additional relation D ~ K z / 2 w 2 , which was verified with the direct computations of 
the Arnold diffusion in systems with a few degrees of freedom. In fact, the value of w s is determined in 
10; 1 1J via the computation of the period of oscillations around a separatrix layer of a nonlinear reso- 
nance that is related to the computation of LE. Due to this indirect method, Chirikov and Vecheslavov 
were able to obtain the variation of the Arnold diffusion constant Da by 50 orders of magnitude! On 
a scale of first 30 orders of magnitude the decay of the diffusion constant D is well described by the 
power law with vr> = 6.5 (see Fig. 1 in |11|). The main message of these amazing calculations is a 
non-exponential decay of D, and hence of the chaos measure w s , with the decrease of nonlinearity 
parameter K . This result is in a drastic difference from the asymptotic Nekhoroshev-like estimates 
based on the KAM theory J6][7] which give exponential decrease of D and w s as K — > 0. Of course, 
there is no formal contradiction since the results for fast Arnold diffusion [11] are always obtained at 
small but finite K values. However, an algebraic decrease with K indicates on an existence of weak 
chaos component with relatively large measure. The heuristic arguments for this phenomenon were 
presented in [11]. According to the results of [11] one has for model A: 

D ~ K^ 2 w 2 s , w s ~ K"' , v a w 2.5 , v D = 2v s + 3/2 , (14) 

for K > 1.6 ■ 10 -5 . Here, w s is a dimensionless measure of the chaotic separatrix layer of the resonance 
between two nearby oscillators. For the range 2-10~ 6 < K < 1.6 TO -5 the decay of D is compatible with 
the power law D cx K 15 but this range of K variation is not very large. The global dependence D(K) 
is fitted by the dependence of Eq. (5.8) in [TT] which however has no complete theoretical explanation. 

The reason why one can hardly compute the diffusion coefficient at smaller K is clear from the 
inspection of the dependence of the variance on time in Fig. [T3] For small K one observes a normal 
diffusion only when the variance exceeds rs 1, below this value the diffusion looks like anomalous 
one with the variance proportional to a power of time. This means that a "random walk" inside the 
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periodicity cell [0, 2ir) is highly correlated, while only cell-to-cell walk demonstrates a normal diffusion. 
For small K the mean first passage time to the next cell becomes extremely large - nearly 10 9 for 
K = 0.03, while for K < 0.03 this mean passage time is of order or larger than the total integration 
time and only the anomalous diffusion is observed. 




Fig. 12 Left panel: Variance <(p(T) - p(0)) 2 ) as a function of time T calculated in a lattice of length L = 64. 
From top to bottom: K = 1.5, 1., 0.5, 0.2, 0.11, 0.08, 0.06, 0.05, 0.04, 0.03, 0.01, 0.005, 0.002. Right panel: 
dependence of the diffusion constant on K in range 0.03 < K < 1.5. Dashed line shows relation D — 2QK 6 ' 5 . 



The obtained properties of diffusion should be contrasted to the properties of LEs, as both quantities 
give some characteristic times of the system. We have demonstrated that these times become extremely 
different for small non-integr abilities, as the Lyapunov exponent A ~ K 1 / 2 decreases rather weakly with 
K while the diffusion constant D ~ K 6 5 drops much more rapidly. We interpret this as indication that 
chaos is mainly "local" , not leading to large deviations of variables. This picture corresponds well to 
the discussed above effective resonances as the origin of chaos: in the triple resonance described above 
in Section [4J the sum of all momenta is a conserved quantity, so that the chaotic dynamics like in 
Fig. [TT] does not lead to a large deviation of momenta involved in the resonance. Indeed, there is strong 
chaotic dynamics inside the triplet resonance, but the sum of three resonant actions is a constant in 
the resonance approximation that would give a zero diffusion coefficient D — 0. However, the resonant 
approximation is not exact and it is destroyed by nonresonant terms and higher order perturbations 
that leads to a finite value of the diffusion D ~ K 6 5 . A mixture of strong chaos, which is however 
bounded due to an additional integral of motion, and a slow but unbounded diffusion produced by 
weak chaos makes the numerical computation of the diffusion rate a rather difficult task. In fact, usual 
very powerful methods discussed in [43], which allowed to compute as small diffusion rate as 10~ 22 , are 
not working in such a situation and only computations at very long times allow to determine directly 
the value of D. 

The physical origins of the power law decay of the diffusion rate with K (fT3"l) are still to be 
understood. The theoretical heuristic arguments presented in |llj assume that in the regime of weak 
chaos a trajectory follows mainly those chaotic resonant layers which have locally most large size. An 
optimization over various resonances leads to a certain power low decay for w s and D which gives 
v s = e = 2.718... and vu — 1.5 + 2e — 6.936... respectively (we remind that for the Chirikov standard 
map w s oc exp(— -k 2 j\j2K) [2"llll)b This theoretical value of the exponent vjj is in a satisfactory 
agreement with the numerical value found at not very small K values. However, at very small values of 
K < 10~ 5 such arguments should be modified to fit an unknown dependence of resonance amplitudes in 
high orders of perturbation theory . According to the heuristic arguments |llj the main contribution 
to diffusion is given by the resonances with an effective resonance harmonic numbers Mq — ln(l/VK) 
with a dimensionless measure of chaos inside one given resonance separatrix layer ws ■ We may argue 
that the number of such layers grows with Mo at least as Mq so that the total measure of weak chaos 
can be estimated as P wc h oc M 2 w s oc (In K) 2 K 2 5 . This dependence is in a satisfactory agreement with 
the data of Fig. \E\ (see the dashed curve there) and the empirical exponent value v wc h ~ 1.6 in ©. 
Thus we can say that our data for the measure of weak chaos are in a satisfactory agreement with the 
numerical results [TT] . 
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On the other hand, the origin of such a weak chaos component is still to be clarified. Indeed, the 
studies and arguments presented in [TTJ did not take into account the strong chaos based on triplet 
resonances which exists at arbitrary small K. This strong chaos component emerges as the result of 
triple primary resonances but it is clear that a similar mechanism can work for higher order resonances 
which may be at the origin of the weak chaos component. On the other hand, the triple-like resonances 
of higher order in K should lead to appearance of a certain number of trajectories with the LEs 
A oc K m / 2 with m > 2 that is, however, is not visible in the distribution of LEs in Figs. I6I7I8I It is 
however, possible that other tiny chaotic layers hide such contributions. Further studies are required to 
clarify these points especially in the regime with large L 3> 15. An indication on the complex internal 
structure of weak chaos provides Fig. [5] above, which demonstrates how a trajectory visits regions with 
different LEs along a very long evolution. 



6 Spreading of chaos 

Above we discussed the local properties of chaos computing the Lyapunov exponents and the diffusion 
rate in the regime when all nonlinear oscillators are populated in the initial state. Another type of 
question appears for the model C2 Q when only one of few nearby oscillators are initially excited with 
the total energy E to t = 1 and 7 = 1 while all other oscillators have zero energy. Since the total energy 
is conserved we face the question on a possibility of energy spreading over the whole lattice of size L. 
This is related to the question of ergodicity of large finite lattices at small energies. In the case when 
both nonlinear terms in the Hamiltonian (the local potential and the coupling) have the same power 
(e.g. the coupling has power 4 instead of 6, such a model can be called model C44) then it is known 
that a thermalization takes place at arbitrary small total energy according to the arguments given in 
pID] . Of course, the time for such global ergodicity grows as a power of system size L. For models with a 
nonlinear destruction of the Anderson localization, we have the terms with powers 2 for local potential 
and 4 for coupling in (J3J, which we call model C24. In this case it is found that a slow subdiffusive 
spreading over the lattice takes place up to very long times t ~ 10 9 (see details in recent papers [TBI 
19,20,24.25,26,27,29,40 ). The model C2 corresponds to a new situation for energy spreading when 
the unperturbed integrable Hamiltonian is nonlinear and the coupling between nonlinear modes has 
higher nonlinearity. In contrast to the FPU problem, here the coupling between modes is local and the 
randomness in local nonlinear frequencies rjk excludes any proximity to a full hidden integrability. 
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Fig. 13 Left panel: Spreading of the second moment (Ak) 2 = ^2 k (k — ko)Ek/ Ylk ^ k anc ^ participation number 
P = Q2 k Eh) 2 1 J]] fc E 2 (inset) vs. time for initial single site excitation fco = in model C2. Data up to t = 10 7 
(empty circles) has been averaged over 1000 realizations of disorder and logarithmic time windows. Long-time 
values until t — 10 8 (full circles) were averaged over 24 realizations. Dashed lines show subdiffusive growth 
(Ak) 2 ~ t a and P ~ t a/2 where the fit of the asymptotic behavior (t > 10 5 ) gave a = 0.55 ±0.01. Right panel: 
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Let us assume that in model C2 with the above local initial conditions the energy spreads over the 
whole lattice of L oscillators with an approximate energy equipartition over L sites. After a rescaling 
of variables of this final state to a new time r — > L 1 / 4 t we come to the model C2 with 7 ~ l/yl and a 
homogeneous initial condition, discussed in the previous sections. In general, the probability of strong 
chaos in such a case scales as P c h oc 7L oc y/~L so that we expect local strong chaos to occur almost 
surely in a sufficiently long lattice. The same it true for the probability of weak chaos even if in this 
case the sum value of the exponents in L is close to zero. Although the probability to observe chaos 
is high, it is important to note that this chaos is mainly local: some modes are chaotic, e.g. triplets 
discussed above, but other modes generally oscillate nearly quasiperiodically. Indeed, in a system with 
many degrees of freedom some modes can be chaotic while others can be close to integrable ones, 
without any contribution to the maximal LE. Thus, it is not obvious if the local strong chaos can allow 
spreading from initial local state over the whole lattice. 

Let us present here simple estimates on the possible rate of such a spreading using results for the 
diffusion in the weak chaos component. We assume that a chaotic spreading populates the number 
of modes N at time t. Using rescaling given above we can argue that the new mode N + 1 will be 
populated due to the weak chaos diffusion after a time scale tg/N 1 / 4 ~ 1/D(j) ~ ^~ v ° ~ N"d/2^ 
This gives us an effective local diffusion rate in N with N 2 /t ~ l/t s ~ \/N( 2u D+ 1 )/ i leading to the 
subdiffusive growth of the second moment TV 2 : 

N 2 ~t a , a = 8/(9 + 2^) . (15) 

For vo = 6.5 we obtain a = 0.3636. However, our results for spreading, shown in Fig. 1131 give 
approximately a — 0.55 that corresponds to vd ~ 2.77. We explain this difference in the following 
way. At the maximum time t max = 10 s , reached in our numerical simulations, the energy spreads over 
a number of modes N ~ tmax so that we have an effective 7 ~ 1/V~N ~ 0.02 which is only at the 
beginning of the decay with the exponent vd shown in Fig. I12f right panel), if we assume a simple 
relation 7 = K , which however still may have an additional numerical factor. It is interesting to note 
that the case with Vd = corresponds to independence of D on iV after rescaling that is the case for 
nonlinear model C44 (with both potentials having power 4 in ([3])) where the spreading goes indeed 
with the exponent a = 8/9 as it is shown in [40] . 




UJ U) u 

Fig. 14 Diffusion D z k(u}) in z-variable for particles k — 1, 8, shown by different color symbols, in model CI 
at L = 8 and 7 = 1.0965 ■ 10" 2 (left), 1.5849 • 10" 3 (middle), 1.2023 • 10 -4 (right). The straight line shows the 
dependence D z k oc exp(— lj/^/j). 

An indirect support to the view point according to which at t max = 10 8 we still did not reach the 
asymptotic spreading exponent a = 0.3636 is based on the numerical computation of the diffusion 
rate in an additional effective degree of freedom described by the equation dzk/dt — sin(wt), where 
qk are dynamical variables in model CI ([3]). Solving these equations in parallel with the dynamical 
equations of motion for q^ we determine the effective diffusion constant D z k(to) for each particle k at 
L = 8. To suppress regular quasiperiodic oscillations we use the window averaging method described in 

[43] computing first the average z k (j) = J^ +1)T z k {t) sm 6 (2Trt/T)dt/ J^ +1)T sm 6 (2Trt/T)dt over time 
interval T = 10 6 and determining the diffusion for each k via the relation D z k(uj) = Y^j>>j>i(zk{j) ~ 
^fe(j')) 2 /((i' — The computation is done for one trajectory with total time t = 10 7 . The initial 

particle energies are chosen to be E k = pf,/2 + q 4 /4: s» 1 at q k = 0. At 7 = we have the particle action 

h = r{l/A)E 3 k /4 /(2^r(7/)) and nonlinear frequency uj k = dH/dI k = S^r^f^E 1 / 4 /(3r(7/8)) « 
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1.2Ej[ w 1.2. The dependence of D z k on frequency u> is shown in Fig. [14] for all k values. In fact, 
D z k(w) gives us the spectral density of an effective noise produced by dynamical chaos. According 
to the results obtained for the modulational diffusion [33], the spectrum of D z k(u) is expected to 
have a plateau of width Aco centered at the resonance u>k ~ 1.2, followed by an exponential drop 
D z k(ui) cx (l/Zkj) exp( — \uj — uik\/Aui). In the picture of triplet resonance we have Au ~ y/j. The 
data of Fig. [14] are in a satisfactory agreement with such a picture showing a decrease of the plateau 
size with the decrease of 7. The plateau is followed by an exponential drop. However, at 7 ~ 0.01 the 
spectral width Auj is still rather large being comparable with ~ 1. At such spectral width even 
the oscillators that are not directly involved in the triplet resonance still will be affected by it. This 
is probably the reason why up to 7 ~ 0.02 we have the spreading of chaos with the exponent a ~ 0.6 
corresponding to a usual diffusion D cx 7 2 in model C2. At 7 ~ 10~ 4 the spectral width becomes 
notably smaller than unity but one needs to go to enormously large times t rnax ~ 10 26 to reach such 
effective values of 7 during spreading of chaos. The value X^7^10~ 5 where there is a change in 
the dependence D{K) detected by Chirikov and Vecheslavov (see Fig. 1 in [11]) would require times 
at least as large as t max ~ 10 33 . Definitely such times remain out of reach of modern computations. 

On the basis of presented results and discussions we can say that the spreading of chaos over the 
nonlinear oscillator lattice of model C2 © goes in a subdiffusive way (|15p with the exponent a ~ 0.55 
up to times t ~ 10 8 . In view of the result of Chirikov and Vecheslavov for the fast Arnold diffusion 
(|13| fTT] it is possible that the exponent will go down toaw 0.36 at times t > 10 26 . The properties of 
chaos spreading behind times 10 33 remain absolutely unknown. During this anomalous slow growth of 
the wave packet size, the chaotic spreading follows the Arnold web of tiny chaotic layers propagating 
mainly along mostly thick ones. However, from time to time a trajectory can go inside thinner layers 
that leads to a strong drop of local diffusion and propagation rates, as well as a significant drop of LE 
(see, e.g., Fig. [§]). It is quite possible that in this regime the energy Ek distribution over the populated 
modes N(t) is still more or less homogeneous, as it is seen in Fig. [JJJ] however, we expect this state to 
be not ergodic within these N modes since chaos is presumably confined inside some "porous medium" 
of Arnold web with very complex structure and topology. In course of spreading, the energy per excited 
oscillator goes down to zero, so that such a process can be considered as an unusual non-ergodic cooling. 

7 Slow diffusion in Hamiltonian systems as deterministic rheology 

The spreading of chaos discussed above goes in a very slow way. In this section we explore a parallel 
with slow rheology processes characterized by small values of the Deborah number [33J 



where t r is a time scale of local relaxation process and t a b s is a time of observation. The values of 
Dr <C 1 correspond to a liquid-like phase while Dr ^> 1 appears for the solid phase. At our initial 
state with one or few excited oscillators we have the relaxation time to be comparable with the inverse 
LE t r ~ 1/A ~ 1, while the observation time of spreading is t bs ~ 10 8 for our numerical simulations. 
Thus we have extremely small values of Dr ~ 10~ 8 for our studies. The parallels with rheology 
processes, which are actively studied in soft matter and porous materials (see e.g [451146] ). can be build 
on the basis of the following arguments: a)in rheology the flow processes are characterized by small 
Dr values that is exactly the case for chaos spreading in model C2; b)often a spreading in a porous 
media is described by a nonlinear diffusion for a density p(x, t) |47j : 



and it was shown recently that this equation gives a good phenomenological description of chaos 
spreading in nonlinear lattices; [26]; c)the Arnold web of chaotic resonance layers forms some kind 
of a porous media along which energy can spreads to larger and larger sizes. Recent experiments 
on gel formed by attractive colloidal hard spheres, suspended in an aqueous solvent, show that the 
spreading of gel is indeed well described by such type of a nonlinear diffusion equation (|17[) with a 
nonlinear flux term [48j . The theoretical models of rheology flow try to explain such a spreading by 
phenomenological statistical models with disorder and metastability (see e.g. [351l50| h In contrast to 
such statistical models, our "rheology" of chaos spreading has purely dynamical and deterministic 
origin. 




1/(AU) « 1 , 



(16) 



dp/dt = D d(p a dp/dx)/dx 



(17) 
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The value of Dr given above should be considered as a global simplified estimate. It is also important 
to see how Dr varies with time t b s of spreading duration. For the model C2 we have A ~ Ik ~ E^ 4 ~ 

TV" 3 / 4 - t^ s a/8 and hence from © we find - l/f^ s 3a/8 oc l/t° fc 7 / > 1. Thus in this model 

D# — > at t bs 00 that argues in a favor of continuation of spreading at infinitely large times. The 
same criterion applied to the DANSE model, which describes the Anderson model with nonlinearity 
/3|V>| 2 and was studied in [IS], gives A ~ I ~ /3/iV - P/t"{* and thus D R - l/(Pt^ a/2) ) still goes 
to zero in the limit of large times (a sa 1/3 for DANSE). The above arguments show that for the 
nonlinearity /3\ip\ 2a studied in [55] we have A ~ I a ~ l/i°°/ 2 and £> fl - l/t^ s " a/2) — > even for 
a = 2,3 (see corresponding values of a given in [26]). Indeed, the numerical results of [26] show an 
infinite spreading for such values of a. 

The above discussion shows that weakly nonintegrable many-body Hamiltonian systems give a new 
interesting example of rheology of chaotic dynamics. These systems are ruled by purely deterministic 
and rather simple Hamiltonian equations of motion. Exploring further statistical properties of such 
a deterministic rheology, generated by Hamiltonian many-body dynamics, is an important task for 
future studies. 



8 Conclusion 

In this paper we studied properties of non-integrable Hamiltonian lattices focusing on the regimes of 
very weak non-integrability. Our main results are scaling relations for the probability to observe strong 
chaos with the largest possible Lyapunov exponent. This probability is proportional to the product of 
the coupling parameter and the lattice length, while the Lyapunov exponent scales as a square root of 
the coupling constant. This behavior is explained by the observation that strong chaos is mainly due 
to resonances that appear when three neighboring sites occasionally have close frequencies. Because 
both the frequency mismatch and the characteristic time scale of the resonance are proportional to a 
square root of the perturbation parameter, the relations above directly follow from this scaling. 

Furthermore, we confirm previous calculations showing that the diffusion time scale at weak non- 
integrability is much larger than the inverse Lyapunov exponent, and relate this to a weak diffusion 
inside the weak chaos component. The measure of this component decreases only algebraically with 
the strength of nonlinear coupling between nonlinear oscillators. The obtained results are in a good 
agreement with the fundamental finding of Chirikov and Vecheslavov [ MTTJIfTT] who first discovered 
this regime, with only algebraic decrease of the measure of chaos and diffusion rate at rather small 
perturbations, and named it the fast Arnold diffusion. 

We also studied the spreading of chaos in such coupled nonlinear lattices showing that the spreading 
goes in an anomalous subdiffusive way. The link between the exponent of this spreading and the fast 
Arnold diffusion are also determined. 

As already mentioned in the introduction, one has to distinguish weakly nonlinear and weakly non- 
integrable systems. There is, however, some analogy between the dynamics of weakly nonintegrable 
lattices studied in this paper and random lattices with weak nonlinearity [2l?ll3l71l2"?U18llil)] . We consider 
homogeneous lattices, where resonances appear randomly due to random choice of initial conditions. In 
random weakly nonlinear lattices resonances are determined by a lattice disorder. So in both cases one 
can expect that chaos is mainly sitting on resonances. For nonlinear homogeneous lattices, resonances 
can "move" as the energies on different lattice sites vary, while in weakly nonlinear disordered lattices 
the resonances are due to disorder and thus are "pinned" . The properties of chaos spreading in the 
latter case require separate investigations. 
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